Exciton transport in thin-film cyanine dye J-aggregates 
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We present a theoretical model for the study of exciton dynamics in J-aggregated monolayers of fluorescent 
dyes. The excitonic evolution is described by a Monte-Carlo wave function approach which allows for a unified 
description of the quantum (ballistic) and classical (diffusive) propagation of an exciton on a lattice in different 
parameter regimes. The transition between the ballistic and diffusive regime is controlled by static and dynamic 
disorder. As an example, the model is applied to three cyanine dye J-aggregates: TC, TDBC, and U3. Each of 
the molecule-specific structure and excitation parameters are estimated using time-dependent density functional 
theory. The exciton diffusion coefficients are calculated and analyzed for different degrees of film disorder and 
are correlated to the physical properties and the structural arrangement of molecules in the aggregates. Further, 
exciton transport is anisotropic and dependent on the initial exciton energy. The upper-bound estimation of the 
exciton diffusion length in the TDBC thin-film J-aggregate is of the order of hundreds of nanometers, which is 
in good qualitative agreement with the diffusion length estimated from experiments. 



I. INTRODUCTION 

In organic materials, excitons, quasiparticles of bound 
electron-hole pairs, act as the intermediates between light 
(photons) and charge (electrons and hole). Understanding 
which physical properties make certain molecular aggregates 
optimal for exciton transfer is one of the main current techno- 
logical goals in organic material research. 

In this article, we develop a computational model and em- 
ploy it to explore excitonic energy transport in a particular 
class of organic materials: cyanine dye J-aggregates. 

Discovered over 50 years ago [1, 2], J-aggregates are 
typically formed by organic fluorescent dye molecules and 
can be identified spectroscopically by the narrowing and ba- 
tochromic shift (J-band) of the lowest electronic excitation rel- 
ative to the monomer band [3, 4]. These structures are charac- 
terized by the unique properties associated with their J-band: 
a large absorption cross-section, short radiative lifetimes, a 
small Stokes shift of the fluorescence line and efficient energy 
transfer within the aggregate [3], which can be used for de- 
signing state-of-the-art photonic devices. 

J-aggregates have been studied both experimentally [5-10] 
and theoretically [11-15] and their applications range from 
being used as "reporter molecules" in mitochondrial mem- 
brane potentials in living cells [16], to photosensitizing sil- 
ver halides in photography [17]. Moreover, J-aggregates are 
imployed in dye-sensitized organic solar cells which provide 
several advantages over inorganic solar cells [18, 19]. Re- 
cently, cyanine dye J-aggregates have been combined with op- 
tical cavities [20-22] or coupled to quantum dots [23, 24] to 
form hybrid systems. 

However, the current understanding of exciton transport 
properties, even for the most ordered J-aggregates, is rather 
limited. This limitation arises from the challenges encoun- 
tered in the experimental characterization of their structure 
[3] and from the lack of information on the dissipation pro- 
cesses. In an effort to overcome these difficulties, the aim of 
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our study is to provide a theoretical model, with a minimal 
number of phenomenological parameters, which is useful for 
the determination of the J-aggregate structure and able to de- 
scribe the exciton dynamics in hybrid excitonic-photonic and 
exci tonic-electronic devices. 

J-aggregates can be found in various structural arrange- 
ments including one-dimensional, planar and cylindrical [3] 
aggregates each exhibiting different optical and exciton en- 
ergy transfer properties [25-27]. The structure of liquid- 
crystal cyanine dyes was initially studied using absorption 
and fluorescence spectroscopy by Scheibe and Kandler [28], 
and more recently using X-Ray diffraction and NMR by Har- 
rison et al. [29, 30]. Nonetheless, the packing structure 
of J-aggregates remains unknown [3]. Different theoretical 
packing models for two-dimensional arrays of these pseudo- 
cyanine dye (PIC) aggregates have been proposed by Naka- 
hara and Kuhn [31]. In this study we focus on modeling 
the exciton transport properties of two-dimensional (2D) thin- 
films using cyanine dye J-aggregates such as those realized ex- 
perimentally in Ref. [32]. These highly efficient light absorb- 
ing thin-films are employed in various opto-electronic systems 
for applications such as lasers and optical switches [33-35]. 
Exciton dynamics in these films, in general, possesses both 
ballistic (quantum) and diffusive (classical) regimes and can 
be analyzed at different levels of approximation. 

Experimentally, transport properties, such as diffusion co- 
efficients of excitons in organic materials can be obtained 
using indirect methods only. These include exciton-exciton 
annihilation [10, 22, 36], photoluminescence quenching [37, 
38], transient grating [39] and photocurrent response [40, 41]. 
In a recent study by Akselrod et. al. [22] the singlet exciton 
diffusion length in a 2D cyanine dye film has been estimated 
to be of the order of 50 nm at room temperature which is more 
than twice larger than the diffusion length measured in stan- 
dard organic semiconductor films [38]. 

Theoretically, exciton transport has been studied using a 
classical hopping model [42]. However, this approach is ap- 
plicable only in the weak Forster coupling regime [43], where 
the exciton mobility is low. Beyond this regime, the tight- 
binding Hamiltonian with classical noise model, proposed in 
the 70's by Haken, Strobl, and Reineker [44, 45] allows for 
a unified description of ballistic and diffusive exciton dynam- 



ics. For perfect structures with translational symmetry this 
model can provide analytic solutions for the moments of the 
exciton wave function [46] that characterize exciton transport. 
Variations of this type of analysis have been provided by oth- 
ers [47-51]. Incorporating a detailed computational ab-initio 
approach into this model increases its accuracy by providing 
more insight on the specific characteristics of exciton trans- 
port and by removing the limitation on the structural symme- 
try. 

To demonstrate the applicablity of our model, we 
present a detailed study of exciton diffusion in three 
types of cyanine-dye J-aggregates, namely TC (5,5'- 
dichloro-3,3'-disulfopropyl thiacyanine), TDBC (5,6- 
dichloro-2[3-[5,6-dichloro-l-ethyl-3-(3-sulfopropyl)- 
2(3H)-benzimidazolidene] - 1 -propenyl] - 1 -ethyl-3-(3- 
sulfopropyl) benzimidazolium hydroxide), and U3 (3-[(2Z)- 
5-chloro-2-[((3E)-3-{ [5-chloro-3-(3-triethylammonium- 
sulfonatopropyl)- 1 , 3 -benzothiazol-3 -ium-2-y 1] methylene } - 
2,5,5-trimethylcyclohex- 1 -en-l-yl))methylene]- 1 ,3- 
benzothiazol-3(2H)-yl] propane- 1 -sulfonate. The devel- 
oped theoretical model provides a tool for the analysis of 
transport properties and can be utilized in modeling hybrid 
excitonic-photonic devices. 

Our findings indicate anisotropy in the exciton diffusion, 
the presence of coherent dynamics at times shorter than tens 
of femtoseconds and finally a dependence of transport on the 
specific molecular excitation parameters. 

The paper is organized as follows. In Section II, we de- 
scribe the theoretical model for the exciton dynamics. In par- 
ticular, we introduce the Hamiltonian of the system and the 
associated Langevin equation. The static and dynamic noises, 
which represent the different types of disorder present in J- 
aggregate films, are discussed. The description of the model 
is completed with a derivation of the diffusion equation. Sec- 
tion III includes details of the calculation of the Hamiltonian's 
parameters and also gives an overview of the Monte-Carlo 
Wave Function method (MCWF) employed in the study of ex- 
citon propagation. In Section IV we analyze exciton transport 
in thin-film J-aggregates of three different cyanine dyes: TC, 
TDBC, and U3 (structures are shown in Fig. 1). We conclude 
the study by summarizing our results in Section V. 



II. THE MODEL 

A. Hamiltonian and single exciton dynamics 

We apply the general exciton theory developed previously 
for molecular crystals and molecular aggregates [52, 53] to a 
specific system - a 2D monolayer J-aggregate of fluorescent 
dye molecules. The exciton Hamiltonian of a 2D aggregate 
can be constructed starting from that of single monomers by 
explicitly incorporating intermolecular couplings. The multi- 
exciton space of the aggregate consists of independent exci- 
ton manifolds that are defined by a specific number of ex- 
citations in the system. The manifolds are coupled to each 
other by exciton relaxation, annihilation or creation processes. 
A more detailed discussion on the single-molecule Hamilto- 
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Figure 1. Structure of the monomer dye molecules TC, TDBC, 
and U3 which form the aggregates. The full lUPAC names of the 
molecules are given in the text. 

nian and the exciton-manifolds is provided in Appendices A 
and B. Here we constrain the exciton dynamics to the single- 
exciton manifold and therefore assume that the exciton density 
of the system is low. The Hamiltonian for a single exciton in 
a molecular aggregate can be written as 

^ _|_ ^cl— bath jjhath 

where °' is the system Hamiltonian which includes the elec- 
tronic degrees of freedom, v^ci-bath system-bath inter- 
action Hamiltonian and H^'^^^ is the bath Hamiltonian. The 
electronic Hamiltonian in the site basis can be expressed as 

N I ^ 

H-'^J2^n\n}{n\ + - J2 Jnm\n){m\, (2) 

n—l n,m—l 

where e„ are the energies of the electronic excitations at each 
site n and J„„j are the couplings between electronic transi- 
tions of monomers at sites n and m, for an aggregate of N 
monomers. The energies e„ are assumed to be equal for all 
sites and disorder in these diagonal energies will be included 
as described in Section IIB. In Eq. 2, \n) = |0...1„...0) cor- 
responds to the state where an exciton is localized at the n-th 
molecule and all other molecules are in their ground electronic 
state. In the aggregate, each monomer is modeled as a two 
level system and the environment is assumed to be a harmonic 
bath formed by the intra and intermolecular vibrations 
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where q runs over all vibrational modes. The system-bath in- 
teraction term in the linear coupling limit is 
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where h\^^ and hnq are the (bosonic) creation and annihilation 
operators for the bath modes at site n and Kg is the coupling 
constant between the q — th vibrational mode and the elec- 
tronic system, assumed to be equal for all sites. 

In general, there are two possible contributions to the cou- 
pling terms: Dexter [54] and Forster interaction [43]. If the 
wave-function overlap between interacting molecules is not 
negligible, the Dexter interaction may give a significant con- 
tribution to the Jnrri Coupling term. While this term is domi- 
nant for triplet exciton transport, singlet exciton transport can 
be accounted for mostly by Forster coupling. Therefore, in 
this case, the major contribution to the J„m coupling terms in 
Eq. 2 is due to the Forster interaction. For a pair of molecules 
m and n, it is given by 
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where (r|(/)„g) and (r|(/)„e) are respectively the electron wave 
functions of the ground and excited states of the n-th molecule 
in coordinate space, Sq is the vacuum permittivity and e is the 
permittivity of the medium. In the following discussion, we 
will assume e = 1. In practice, numerical calculations of 
the Forster term using Eq. 5 can be computationally heavy, 
especially for large structures. For J-aggregates, where the 
distance between stacked molecules is comparable to the spa- 
tial extent of each molecule, it is possible to use the extended 
dipole model [55]. Within this model, the Forster interaction 
between two particular electronic transitions is parametrized 
by a transition charge, q, and a transition dipole length I. The 
interaction term can therefore be simplified and written as a 
sum of Coulomb interactions between the transition charges 
located on different molecules. 
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where r+,^ is the distance between the charge +q located 
on the n-th molecule and the charge —q located on the m-th 
molecule. 

Fluorescent dyes may self-aggregate in a number of dif- 
ferent structures [56]. When dealing with two dimensional 
films of cyanine dye aggregates, the brickstone model is one 
of most commonly employed models which can account for 
the experimentally observed optical properties of these ag- 
gregates. Therefore, in our work, the molecular arrangement 
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Figure 2. Brickstone lattice of a 2D J-aggregate. Each brick repre- 
sents a molecule. Green arrows indicate extended molecular transi- 
tion dipoles and I their length. The three lattice parameters are 9, Lx 
and Ly which can be identified by the lattice vectors and Ly. 



within the 2D layer is modeled as a brickstone lattice [7], as 
shown schematically in Fig. 2. 

In this model, the dye molecules are stacked parallel to 
each other and subsequent rows are displaced by an angle 9. 
The lattice is characterized by two lattice vectors and hy. 
Given the Hamiltonian of the system, we proceed to derive 
an equation for the quantum evolution of the system's wave 
function on this lattice. 

The quantum dynamics of excitons can be numerically sim- 
ulated using various methods, including density matrix evolu- 
tion schemes [57, 58] and also diffusion [59, 60] and quantum 
jump [61, 62] wavefunction approaches. The latter are based 
on propagating the system's wavefunction rather than the den- 
sity matrix and are therefore computationally more convenient 
to be employed for modeling large systems. 

We define a quantum stochastic equation following the 
Langevin procedure, as described in Ref [63]. Before includ- 
ing any source of noise or relaxation, the equation of motion 
for a single exciton in the aggregate is simply the Schrodinger 
equation 



d\m) ^ ^fj. 

at h 
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where \i'{t)) is an exciton wave function in the single exciton 
manifold and is defined in Eq. 2. The initially excited 
exciton state will eventually decay back to the ground state 
due to interaction with the environment. Processes such as 
population relaxation and dephasing between excited states 
and the ground state can be included as follows 
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dt 
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where the 



Hamiltonian H'^^'^ 



cl 



effective 

T fi (pinj includes a decay term and a stochastic 
fluctuation term ry^ (t) for site m and channel /i which repre- 
sents a dynamic noise force and is introduced to conserve the 
norm of the wave function. The (7^ are Lindblad operators 
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for two channels /i € {l,(f>} (relaxation or dephasing) for 
each site m. In our case, these terms can be expressed as 



(9) 



where describes single exciton relaxation and Cf^^ is in- 
troduced for dephasing processes, which are associated with 
the relaxation and dephasing F*^ rates. The states \m) span 
over the single exciton manifold, and the state |0) denotes the 
ground state, where no exciton is present. The relaxation and 
dephasing rates are taken to be equal for all sites, because the 
aggregate is constructed of identical molecules which are as- 
sumed to be in identical local environments. It is easy to see 
that with the noise operators written in the site basis this model 
is equivalent to the Haken, Strobl and Reineker model [44]. 



B. Static and dynamic noise 

The environmental noise of the open quantum system needs 
to be introduced into this J-aggregate model in its various 
manifestations. These include fluctuations of site energies and 
of the site-to-site couplings as well as a term which induces 
exciton relaxation. When studying transport properties within 
the single exciton manifold, one is mostly interested in the dy- 
namics occurring on a timescale sufficiently shorter than the 
exciton relaxation time, this corresponds to saying that over 
the time scale of the dynamics, in Eq. 10, Cl^ ~ 0. Thus, the 
exciton number can be assumed to be constant. Moreover, 
for the sake of simplicity in the provided examples we do 
not consider fluctuations of intermolecular interactions. Al- 
though these fluctuations can give a noticeable contribution to 
the transport properties [47] the main effect of noise can be 
captured by the fluctuations of the sole site energies. 

Fluctuations of site energies can originate from low- 
frequency intramolecular vibrations, molecular conforma- 
tions, binding and un-binding events, or charge fluctuations in 
a substrate or solvent which couples locally to the electronic 
states of molecules. Such fluctuations have been successfully 
modeled in natural molecular aggregates using molecular dy- 
namics methods [64, 65]. Detailed studies of the physical ori- 
gin of this noise is beyond the scope of the current paper 
Here, we only emphasize that one can introduce a distinc- 
tion between static and dynamic noise based on the correlation 
time characterizing the fluctuations as compared to the exci- 
ton propagation time. For instance if the conformation of a 
molecule changes on a nanosecond timescale, the fluctuation 
of the site energy remains correlated during the lifetime of an 
exciton which is of the order of tens of picoseconds. There- 
fore, each exciton sees different, yet static local fluctuations. 
In contrast, dynamic noise requires instantaneous fluctuations 
of the site energy with a shorter correlation time, that is in a 
timescale smaller or comparable to that of the energy transfer 
process. 



Static noise can be accounted for in the Hamiltonian, Eq. 
2, by introducing random shifts in the site energies. In our 
model, we extract the random shifts from a Gaussian distribu- 
tion 
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where a is the variance, and eo is the transition energy of each 
isolated molecule. The static noise distribution is assumed to 
be identical for all monomers. Other choices of static disor- 
der models are possible. For instance, a more general Levy 
distribution can be used, which results in different state distri- 
butions and optical properties [66]. 

Dynamic noise directly enters the Schrodinger-Langevin 
equation, Eq. 8 described in the previous section. The re- 
quirement that the norm of | V'(i)) is preserved imposes certain 
properties on the dynamic noise force rjit) [63]. First, the av- 
erage of j]{t) should vanish, i.e., {r]{t)) = 0. Otherwise, one 
has to renormalize the unperturbed Hamiltonian in Eq. 2 to in- 
clude the noise-induced energy shifts. Second, the stochastic 
force ri{t), which is a result of multiple uncorrected micro- 
scopic movements of the molecules constituting the aggregate 
lattice, does not have memory, i.e., {rj{t)rj{t' )) oc V'^5{t — t'). 
This corresponds to the Markov approximation, in the sense 
that the system-bath interaction is assumed to be quasi instan- 
taneous and successive interaction events are not correlated. 
In general, the Markov approximation holds, as long as the 
bath correlation time is much smaller than the time over which 
one extracts properties of the system. Therefore it is only 
necessary that the bath correlator be sharply peaked [63, 67]. 
Within this argumentation the dynamic noise can be character- 
ized by a single parameter F"^, which describes both the bath 
dynamics and the system-bath coupling. While avoiding a de- 
tailed description of the specific environment of a J-aggregate 
we can make the rather qualitative assumption that F"^ is of 
the order of ksT. Such assumption is not strictly justified but 
can be intuitively explained as following. 

If we assume that the bath modes can be modeled as a set 
of harmonic oscillators, we can define, within linear response 
theory, a system bath interaction term linear in the displace- 
ment of the bath modes. It follows that the bath correlator can 
be expressed as [68] 



(11) 

r /h,,,R\ 

dwj (uj) 



coth ^ cos (cjt) — isin (lut) 

where the bath operators are given in the interaction picture 
and we have assumed that the noise correlator is the same for 
all sites and that each bath is uncorrelated from that of other 
sites. The second line corresponds to the approximation of a 
continuous bath spectrum with a spectral density J(a;). Such 
approximation provides a qualitative correspondence between 
the temperature and the dynamic noise in the system. In gen- 
eral, a more quantitative analysis should include explicitly the 
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bath vibrational modes. By using the continuous limit of the 
spectral density we intentionally simplified the model mak- 
ing it independent of the specific molecular vibrational modes. 
Several forms of bath spectral densities are used for modeling 
dissipative quantum dynamics in molecular aggregates [69- 
72]. We choose to employ the Ohmic exponentially cutoff 



spectral density J{uj) 



where A is the reorga- 



nization energy and ujc is the cutoff frequency. This spectral 
density can be integrated analytically in Eq. 1 1 and it also 
has been used for the simulation of exciton dynamics in natu- 
ral molecular aggregates in photosynthetic systems [72, 73]. 
In our model, the correlator in Eq. 1 1 can be used only 
in the limit when the time dependent characteristics of exci- 
ton dynamics remain steady on timescales sufficiently longer 
than the bath correlation time. In this case the temperature- 
dependent dephasing rate can be defined as [72, 74] 



r"* = 27r 



ksT A 
h hujr 



(12) 



where physical properties of the bath and the system-bath cou- 
pling are introduced through the slope of the spectral density 
at zero frequency |^^o = [71]. Eq. 12 is strictly 
vaUd for T > ^ [75]. If we use, for example, values of 
LUc = 150 cm^^ and A = 35 cm^^, typical for the analy- 
sis of quantum dynamics in photosynthetic systems [72, 76], 
the dephasing rate w lAksT. The reorganization en- 
ergy of J-aggregates is comparable to this value, for instance, 
Atdbc = 29cm-i [77]. By setting T"^ = hgT we choose a 
lower bound for the exciton dephasing rates. For room tem- 
perature we thus have F"^ « 26 meV. In reality, there may 
be more sources of dissipation, and different estimates for the 
bath spectral density [71, 73] can give values of the dephasing 
rate, F-^ which are several times larger than the value we use. 



This line of reasoning was first made by Scher and Lax 
[78, 79] for the calculation of the conductivity in doped semi- 
conductors. It was later applied to the calculation of the photo- 
conductivity of organic molecular crystals [48-50] and to the 
bond-precolation model [80]. The method developed by Scher 
and Lax is based on linear response theory [81]. To apply it 
to the excitonic system, we will have to imagine an external 
perturbation which when applied, drives the motion of the ex- 
citons. This can be physically achieved, e.g. through some 
non-resonant coupling with an external field which creates a 
spatially inhomogeneous change in the local potential term in 
the Hamiltonian (cf Eq. 2). However, the actual implemen- 
tation of this scheme is technically challenging if not impos- 
sible. Furthermore, the formalism of Scher and Lax requires 
the density matrix of the excitonic system to be an identity 
in thermal equilibrium; as pointed out in Ref [49], this may 
potentially lead to internal inconsistency in the mathematical 
derivation. To avoid these complications, in the following, we 
will present a simpler way to obtain the relationship Eq. 1 3 
for excitonic systems, without invoking linear response the- 
ory. This approach is essentially the quantum extension of the 
classical approach described in Ref. [82]. 

To get started, we consider a dilute system of excitons 
where their reciprocal interaction can be ignored. We will 
trace the motion of a "tagged" exciton. The probability 
P(r, t) of finding that exciton at location r and time t is given 
by P{r,t) — Tr[|r) (r|p(i)], where p{t) is the density matrix 
of the total system (i.e. exciton plus the bath) at time t. Using 
the following identity 
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C. Diffusion model 

In this section, we consider methods for calculating the dif- 
fusion constant D from the transport properties of the exci- 
tonic system. For classical Brownian motion, it is well-known 
that the diffusion constant is related to the long-time limit of 
the second moment ((r(t) — ro)^) evaluated for the initial 
condition where the particle is localized at a single point Tq 
in space. Explicitly, it is given by 



D 



(13) 



where d is the dimension of the space. For excitonic sys- 
tems, in the first place, we will not assume that exciton motion 
can be described as the motion of classical Brownian parti- 
cles. Here our goal is to show that the relationship above (Eq. 
13) holds even for excitonic systems where a fully quantum 
mechanical treatment is assumed. In other words, the non- 
equilibrium dynamics of an initially localized exciton can in- 
deed tell us about the diffusion constant which describes the 
excitionic motion in (or near) equilibrium. 



where d is the dimension, and f = J^r' r'|r')(r'| is the posi- 
tion operator for the exciton, we write the probability of find- 
ing the exciton as 

oo 

/ dke^'^Pkit), (15) 



where Pfe(t) = Tr[e-''^'-^(*)/5(0)], r{t) = U^t)rU{t), where 



U{t) 



-iHt 



and H is defined in Eq. 1 . 



Since we are considering the fluctuations of the exciton in 
some steady-state (long-time t — > oo) limit, measurable phys- 
ical quantities, including the diffusion constant, should not de- 
pend on the initial condition. We can therefore choose the fol- 
lowing initial state /5(0) =15(8) /5B/Tr(l5), where I5 is the 
identity matrix for the system, i.e. the exciton, and ps — 



/Tr(e 



) is the density matrix of the bath. 



which is assumed to be in thermal equilibrium. Now, inserting 
the resolution of the identity, I5 = |ro) (rg |, we can write 
h{t) - (l/tr(is))Ero e-'''"-''Tr[e-'k-Af(t)|i.o)(ro|®pB], 
where Af(i) = f(i) — Tq. Performing the cumulant expansion 
[63] and keeping terms up to the second-order, we obtain 
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tr(Is) 



ro 



where {[k- Ar{t)]-^)o = Tr{[k- Af(t)]2|ro)(ro|(^PB}. Then, 
we arrive at the diffusion equation 




dephasing quantum jump 



'0 free evolution e~^"''' ' 



*(t + 5t) ■ 



^P{r,t)^ D,,{t)^^P{r,t), ill) 



dt 



where 



^ ^|(Af,(t)Af,(t)) ^ lj^M^J\t) (18) 



is the time-dependent tensor of the diffusion coefficients. 
Here we have defined the second moments M^j ' {t) = 
{A'ri{t)Afj{t)). In order for the diffusion coefficients Dij{t) 
to converge in the long time limit (t — > oo), the right-hand 
side of Eq. 1 8 should scale at most linearly in t. This re- 
sult coincides with that in Eq. 13, which is a special case of 
isotropic diffusion. 



III. COMPUTATIONAL DETAILS 



A. Monomer properties 



Figure 3. Monte Carlo Wavefunction stochastic dephasing jump evo- 
lution. Each stochastic trajectory will evolve either by making 
a dephasing quantum jump or according to a non-Hermitian effec- 
tive hamiltonian . The choice of making a dephasing jump is 
determined by a Monte-Carlo type algorithm, where one extracts a 
random number e and compares it to the dephasing jump probability 
The procedure is repeated up to the final propagation time. One 
then repeats the procedure to obtain an ensemble of trajectories and 
average over this ensemble to obtain the desired observables. 



from IV'(O) according to the stochastic scheme depicted in 
Fig. 3. In particular, at each time step, if the stochastic vari- 
able extracted, e, is larger than the quantum dephasing jump 
probability Q,^, the wavefunction will propagate freely under 
the non-Hermitian effective Hamiltonian introduced in 
Eq. 8 



m + 5t)) 



(19) 



Density-functional calculations of the molecular structure 
and electronic excitation spectra were performed with the 
quantum-chemistry package Turbomole, version 5.10. [83]. 
Triple-C valence-polarization basis sets (def2-TZVP [84]) 
were used together with the hybrid functional of Perdew, 
Burke, and Ernzerhof (PBEO) [85]. Dielectric properties of 
the medium where accounted for by using COSMO [86] as 
implemented in Turbomole. 

To calculate the extended dipole parameters of the fluores- 
cent dyes, the HOMO and LUMO orbitals of the molecules 
were computed on a homogeneous spatial grid. The grid steps 
are dx = 0.5 A, = 0.5 A in the plane of the molecu- 
lar backbone, and dy — 0.25 A in the direction orthogonal 
to the backbone. The Forster, Eq. 5 and the Dexter interac- 
tions between pairs of molecules were calculated for center- 
to-center displacements scanned over x = [—60, 60] A and 
y = [—10, 10] A with step size Aa; = Ay = 0.5 A . 



B. Monte Carlo wave function propagation 

Eq. 8 is a Markovian stochastic open quantum system equa- 
tion, which is equivalent to a Lindblad master equation for the 
density matiix [63]. As such one can evolve single stochastic 
eigenfunction trajectories instead of the full density matrix, 
using the Markov Monte-Carlo wave function method [61]. 

The initial wavefunction [■'/'(O)) = |00...1....0) is defined as 
an exciton localized at the center of the lattice (i.e., in the po- 
sition ; ^] ) . At time t + St, \-)jj{t + St)) can be obtained 



where Af = ^y {ip(t + St)\'>p(t + St)) is the normalization 
constant. However, if e is smaller than il^, a quantum de- 
phasing jump will occur The quantum dephasing jump, a 
specific type of quantum jump, is described as a flip of the 
sign of wavefunction's coefficient corresponding to a site m 
and corresponds to applying the dephasing jump operator C^, 
Eq. 9, to the wavefunction [61]. The phase jump occurs 
in position N — round{-^nxny) where il0 is defined as 

= ^T't'Ti^ny, and where T"^ was given in Eq. 1 2 while St 
is the time step which is assumed to be small enough so that 
ilcf, ^ 1. For the trajectories controlled by exciton dephasing 
only, the effective Hamiltonian is the same as Hq, thus 
the norm of the wavefuction in Eq. 19 is conserved. 

To analyze exciton transport properties of TC, TDBC and 
U3 aggregates the exciton wavefunction was propagated on a 
lattice of 2601 monomers (ux = Uy = 51 ) for a total time 
of t = 100 fs. The quantum trajectories were averaged over 
1000 different realizations of static disorder (convergence on 
the populations was reached with 1000 realizations within an 
error < 4%). The time step in the propagation was set to 
St = 0.6 as thus the probability of the quantum jump within 
the step is sufficiently smaller than 1, and the wave function 
was collected at each femtosecond. The long range interaction 
between the molecules has been accounted for within the cut- 
off distance ^cutoff = 6 • L^, more details about this cutoff can 
be found in Section IV A. The excitation was injected at the 
energy of the maximum of the J-band and at time zero, only 
the central site (26,26) of the lattice was excited. The results 
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were obtained in general over the range F** = [20 — 110] meV 
and (7 = [0—110] meV. Exciton diffusion coefficients where 
estimated from a linear fit of the second moments of the exci- 

(2) 

ton distribution functions as per Eq. 18. 



IV. RESULTS AND DISCUSSION 
A. Model Parameters 

1. Monomer calculations and absorption spectra 

The structures of the TC, TDBC, and U3 cyanine dye 
molecules have been optimized using density functional the- 
ory (DFT). In the computation, we considered single-charged 
anions and assumed that the Na+ ions were dissociated. For 
each molecule, although the conjugated part of the structures 
remain almost planar, there are many conformations which 
differ slightly by the orientations of the sulphonated group 
side chains and are all closely spaced in energy. Examples 
of such conformations for the cK-isomers are shown in Fig. 1. 
The difference between the ground-state energies of cis and 
frani-isomers of the molecules is of the order of hundreds of 
meV. We have chosen cis-isomers as our reference structure 
because it is most likely that in this conformation, with the 
sulphonate groups pointing towards the surface and binding 
chemically or physically to it, that one would obtain the ob- 
served 2D monolayers [32]. For each optimized molecular 
structure, we computed the 100 lowest electronic excitations 
which fall in the energy range 2-7 eV. The computed spectra 
of the molecules shown in Fig. 1 are provided in Fig. 4. The 
strong lowest excitation can be accounted for by the HOMO 
to LUMO transition for more than 98%. Such transition is 
generally assigned to the lowest absorption peak observed in 
the monomer spectra of fluorescent dyes. The second elec- 
tronic transition is separated from the lowest one by about 0.7 
eV, 1.0 eV, and 1.4 eV for TC, TDBC and U3 respectively. 
Moreover, the oscillator strengths of these subsequent transi- 
tions are about two orders of magnitude smaller than that of 
the lowest one. These results support the two-level model pro- 
vided that the static energy disorder is of the order of 100 meV. 
The computed frequencies of the lowest electronic transitions 
are systematically blue-shifted by about several hundreds of 
meV as compared to the experimental values. This shift typ- 
ically occurs in DFT calculations with the PBEO functional 
[87]. Similarly to what was found for the ground state ener- 
gies, the lowest electron transition frequencies for trans and 
c«-isomers differ by hundreds of meV's. 



2. Couplings 

The extended dipole parameters for the lowest electron ex- 
citation are calculated within the Frontier Orbital Approxima- 
tion (FAO) [88], in which one assumes that only the HOMO- 
LUMO transition is involved. To obtain the I and q parame- 
ters in the extended dipole coupling formula Eq. 6, we assume 
that the only type of interaction involved is Forster interaction 
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Figure 4. Computed spectra of electronic excitations in TC, TDBC, 
and U3 c;.s-isomers. The lines are broadened with Lorentzians of 20 
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Figure 5. 2D map of the Forster interaction between two TC 
molecules aligned to each other. The axes show center-to-center dis- 
placement of the molecules. The inset figure shows the orientation of 
the molecules. The color map represents asinh( J^) to emphasize the 
long-range behaviour of the interaction. The blue-green colors define 
the negative frequency shift of the lowest electronic excitation, while 
the yellow-red colors indicate the positive shift. 



(Eq. 5) and fit the interaction between two molecules on the 
X ~ y plane to the Forster results. An example of the calcu- 
lated interaction contour plot for a pair of TDBC molecules is 
shown in Fig. 5. For intermolecular distances larger than 2 A, 
the profile reproduces the interaction of two dipoles and can 
easily be fitted with the extended dipole formula. The largest 
positive shift of the electron transition is obtained when the 
molecules are displaced along the y-axis (direction orthogo- 
nal to the backbone of the molecules). The largest negative 
shift of the electron transition corresponds to the case when 
molecules are displaced along the x-axis approximately by 
half of the length of the molecule. These results are consistent 
with the extended dipole model. The computed propeities of 
TC, TDBC, and U3 monomers are summarized in Table I. We 
also computed Dexter couplings and these were much smaller 
than the corresponding Forster terms for distances of the order 
of the physical spacing between molecules and were therfore 
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Dye n[sY] fi[D] ^tHL [D] I [A] q [e] 

TC 3.3 8.9 8.3 9.1 0.20 

TDBC 2.9 13.1 11.0 10.5 0.22 

U3 2.4 14.6 12.5 11.1 0.24 



Table I. Computed excitation properties of fluorescent dye 
molecules. Q is the frequency of the electronic excitation, fi is 
the transition dipole associated with the transition computed using 
TDDFT, fiHL is the transition dipole computed using HOMO-LUMO 
orbitals only, I is the length of the extended dipole, and q is the charge 
associated with the extended dipole. 



neglected. 



3. Lattice Parameters and absorption spectra 

The brickstone lattice parameters where determined as fol- 
lowing. The horizontal distance between monomers = 
[Lj-l was chosen to be the optimized length of the molecule 
in the c«-geometry /{j*** plus twice the Van der Waals radius 
of the chlorine atom. In particular we chose the longitudi- 
nal Van der Waals radius determined for a C-Cl type bond 
rci = 1.58 A, as reported in Table 11 of [89]. The angle be- 
tween monomers was set to 6 = tan^^ ■ Having fixed 

6 and we then determined the vertical distance between 
layers of monomers Ly -smd by fitting the theoretical position 
of the J-band in the absorption spectra to the experimental re- 
sult. All of these parameters are reported in Table 11 



Dye /ff"[A] L^[k] LysiriefA] 6* [rad] 

TC 15.01 18.17 3.815 0.421 

TDBC 17.36 20.52 4.600 0.404 

U3 19.72 22.88 5.120 0.427 



three dyes. The oscillator strength of a particular transition 
is proportional to the square of the corresponding transition 
dipole. To account for the static disorder the transition fre- 
quencies of the monomers were taken from a Gaussian dis- 
tribution of width 70meV. In Fig. 6 the calculated spec- 
tra of J-aggregates are shown as compared to the energies of 
single molecule excitations. The peak positions where fit- 
ted to the experimental results. To obtain the correct shift 
within 5% of error, we determined that the cut-off distance 
for the molecule-molecule interaction should be not smaller 
than ^cutoff = 6 • L-c. We observe the typical band narrowing 
of the J-band whereas the exact vibrational structure of the 
monomer and the J-band cannot be captured with this simple 
analysis and is beyond the scope of the present paper More- 
over, the present analysis doesn't account for the line shift due 
to the non-resonant (van der Waals) interactions. 
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Figure 6. Calculated electronic excitation spectra of the three molec- 
ular aggregates, TC, TDBC and U3, as compared to the electronic 
transitions of monomers that form the aggregates. Zero frequency 
correspond to the energy of a monomer excitation. The transition 
frequencies of the monomers were taken from a Gaussian distribu- 
tion of the width 70 meV to account for the static disorder. 



Table II. Lattice parameters for the three molecules. Zq" is obtained 
from the DPT optimization of the molecular structure, Ly parameters 
were obtained from fitting the theoretical spectrum to the experimen- 
tal I-band shift. 

The measured energies of the lowest electronic transitions 
for TC, TDBC and U3 dyes in solution and also in the aggre- 
gated form are collected in Table III. 



Dye 


Monomer transition [cV] 


I-band [eV] A [eV] Ref. 


TC 


2.900 


2.613 0.287 [90] 


TDBC 


2.396 


2.115 0.281 [91] 


U3 


1.864 


1.571 0.293 [92] 



Table III. Experimental data for electron excitations in monomer and 
J-aggregated dyes as well as A the shift between monomer and J- 
band. 

To estimate the shifts of the J-band due to the molecular 
aggregation the excitonic spectra of the aggregates have been 
computed by diagonalizing the Hamiltonian, Eq. 2, for all 



B. Quantum exciton dynamics 

The developed model accounts for both coherent and inco- 
herent properties of exciton dynamics that can be monitored 
by following the spatial distribution of the exciton population 



p^At) = {^j\m){m\ij), (20) 

where is the pair of cartesian indices for a particular 

molecule on the 2D lattice and the bar above the expression 
coiTesponds to the ensemble average over quantum trajecto- 
ries. In addition, to characterize specifically the role of coher- 
ences in the exciton transport we analyze the two-point one- 
time correlations between the central site (the point of the ex- 
citon injection) and the current site 

Q,it) = mm)mm\. (21) 

A similar quantity has been used previously in Refs. [93, 94] 
to estimate the exciton delocalization length in natural molec- 
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ular aggregates. If the exciton transport is completely inco- 
herent and represented by the hopping of exciton population 
between sites, the correlation function, Eq. 21, should remain 
zero for all times provided that no initial site-site correlations 
were created. This is consistent with the conventional Bloch 
equations, where the coherence dynamics and the population 
dynamics are separated. In Figs. 7 and 8 we show an ex- 
ample of the population and coherence dynamics in TDBC 
J-aggregate when the dynamic disorder is F"^ — 30 meV and 
static disorder is <t = 70 meV. 

Based on the discussion in Section 11 B, this value of dy- 
namic disorder represents a lower bound to the exciton de- 
phasing rate (upper bound to the exciton diffusion length) 
at room temperature. One can see that exciton transport is 
anisotropic and the population spreads in time following ap- 
proximately an elliptic shape with major axis x and minor 
axis y (such directions are indicated in Fig. 2). The popu- 
lation spreads about 2 to 3 times faster in the x direction than 
it does in the y direction. This can be explained by analizing 
the direction of maximum coupling between nearest neigh- 
bors. In fact, looking at the lattice in Fig. 2 and considering 
the molecule placed at the origin of the purple lattice vectors, 
we see that it has four first nearest neighbors, one of which 
is indicated by the Ly vector The sum of the nearest neigh- 
bors coupling vectors taken in pairs along the x axis is greater 
than that along the y axis. The fact that the exciton is trans- 
ported in the direction of maximum coupling is also seen in 
Fig. 9 where the population is more rapidly transferred to site 
(26,27) along the principle x direction respect to site (24,27) 
which is along the y axis. 

By comparing Fig. 7 to Fig. 8, we notice that coherences 
are spread and suppressed much more rapidly than popula- 
tions. However, the principal transport axis remains the same, 
i.e it is the major axis of the ellipse. This is true for all 
examined values of static and dynamic disorder. At times 
shorter than about 30 fs the dynamics is coherent and inter- 
esting shapes such as a four leaved clover at 2 fs appear in the 
coherence plot. At longer timescales incoherent diffusion pre- 
vails and these features mostly disappear. The peaks of these 
correspond to beats in coherences between second to fourth 
nearest neighbors but such beats die off rapidly over a couple 
of tens of femtoseconds. While quantitative differences are 
observed between the molecules, the general trend is qualita- 
tively similar to the one shown in the example. As a possible 
extension of this initial study of coherent dynamics it would 
be interesting to explore the phase directed exciton transport 
in these two-dimensional aggregates as has been done for the 
one dimensional case [95]. 

While the initial condition (a localized exciton) determines 
a rather large population of the central site even at relatively 
long times, this does not globally affect the exciton diffusion 
properties of the aggregate. Furthermore, the localized exci- 
tation can be viewed as exciton injection from a donor 



C. Diffusion coefficient and diffusion lengtli 

The computed wave function provides the most complete 
source of information about the exciton dynamics. Its second 
moment can capture the main diffusive and ballistic features 
of the transport. For a homogeneous system with stochastic 
dephasing noise andtranslational invariance, the moments of 
the wave function in Eq. 8 can be written analytically [46]. 
If the exciton is initially localized on a particular molecule 
one should expect a ballistic exciton propagation (the second 
moment scales quadratically in time, M'^' cx t^) followed 
by the diffusive motion (the second moment is linear in time, 
Af '^^^ oc t). For systems with static disorder similar transport 
regimes should be observed provided that the dynamic noise 
is strong enough to overcome exciton localization. To verify 
this the second moments of the wave functions each trajec- 
tory were computed over an interval of 100 fs with a time step 
of 1 fs and averaged over a thousand trajectories. In Fig. 10 
an example of Af(^)(i) for the TDBC J-aggregate with static 
disorder a = 70 meV and dynamic disorder F"^ — 30 meV 
is shown. On timescales shorter than about 30 fs the scaling 
of Af*^^Hs approximately quadratic, while the longer time dy- 
namics reflects the diffusive transport. A similar tendency is 
observed for other molecules within the whole studied spec- 
trum of (T and F"^ with the free exciton propagation time 
shrinking down to below 20 fs for large values of F"^ . 

The exciton transport properties depend on the initial in- 
jection energy Ein- Results showing the second moment 
as a function of Ein together with the exciton density of 
states (DOS) are reported in Fig. 11. We see that there is 
a maximum in the moments somewhere between the J-band 
(Ein — —281 meV) and the monomer band (Ein = meV). 
This is true at all times and the position of the maximum is dif- 
ferent for the XX component than it is for the yy component. 
This can be correlated to the large density of states in that en- 
ergy interval (Fig. 1 1 top panel). The high density of states 
is not sufficient to explain the exact position of the maximum, 
in fact, while the density of states has a maximum at about 
— 72 meV, the second moment is peaked at about —170 meV 
for the XX component. The same trend in the second moments 
is observed for TC, and the only apparent difference is that the 
values of the second moment are smaller One aspect to keep 
into consideration is that although there is a large density of 
states around 60 — 80meV, these states mostly have small 
oscillator strengths (See Fig. 6). To investigate further the 
origin of the differences in these maxima one should look into 
the structure and spatial overlap of the exciton states. 

To estimate the exciton diffusion coefficients the linear fit 
was taken over the time interval Ati = 30 — 100 fs to exclude 
the initial ballistic propagation of the exciton. For most val- 
ues of static and dynamic disorder, the boundary effects asso- 
ciated with the finite size of the simulated lattice is negligible 
on the timescale of 100 fs. However, for weak dynamic noise 
the exciton wavefunction reaches the boundary of the simu- 
lated lattice along the y-axis at about 70 fs. In such cases the 
yy component of the diffusion coefficient Dyy was fitted on 
the time range 30 — 70 fs. With the initial condition of injec- 
tion into the J-band, the diffusion coefficients were calculated 



10 




50 60 70 80 90 100 50 60 70 80 90 100 

X coordinate [nm] 



Figure 7. Contour plots of logarithm of exciton population in TDBC J-aggregate projected onto coordinate space at 2, 10, 40 and 80 is. Here, 
the static disorder is ct = 70 meV and the dynamic disorder is V*" = 30 meV. Population spreads more rapidly in the x direction than in the 
y direction. This behavior is observed for all studied values of dynamic and static disorder and can be explained by the couplings between 
monomers as described in the text. 
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Figure 8. Contour plots of logarithm of exciton population in TDBC J-aggregate projected onto coordinate space at 2, 10, 40 and 80 fs. Here, 
the static disorder is cr = 70 meV and the dynamic disorder is F'* = 30 meV. Coherences spread more rapidly in the x direction than in the 
y direction. This behavior is observed for all studied values of dynamic and static disorder and can be explained by the couplings between 
monomers as described in the text. Further, coherences spread and decay much more rapidly than populations, as can be seen by comparing to 
Fig. 7. 



using 

A'4p{t)^2Dut (22) 

for each type of aggregate. The results are shown in Fig. 12. 
Two distinct characteristics emerge from this plot. First of all, 
we notice that diffusion is greater for U3 than it is for TDBC 
which in turn is greater than that of the TC aggregate. This can 
be explained by looking at the physical characteristics of the 
molecules (Tab. 1). U3 has the largest transition dipole, this 
leads to stronger coupling between the monomers and thus to 



more rapid exciton density transfer. The diffusion coefficients 
normalized with respect to the square of the corresponding 
transition dipole (not shown in the figures) are similar for two 
molecules, TDBC and U3, while the normalized diffusion for 
TC is still higher. We attribute this difference to the closer 
packing of TC molecules. This trend is robust to both static 
and dynamic disorder, so long as both are finite and not so 
large as to lead to localization. This implies that indepen- 
dently of temperature, within the studied interval, the trans- 
port efficiency is dictated by the specific physical properties 
of the molecules. Of course, the geometric arrangement is a 
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Figure 9. Populations of sites neighboring source (26,26) as a func- 
tion of time for TDBC with a = 70 meV and T"^ = 30 meV. Popu- 
lation is transferred most rapidly from the central site (26,26) to site 
(25,27), the first nearest neighbor and then to site (26,27) which is in 
the direction of maximum transport. 
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Figure 10. Second moment of the exciton wavefunction in time for 
TDBC aggregate with static disorder cr = 70 meV and dynamic dis- 
order r*^ = 30 meV. Both XX and yy components are included with 
their respective error bars. A transition from a ballistic regime to a 
diffusive regime is observed at about 30 fs. For all studied values of 
static and dynamic disorder these two regimes are observed and the 
transition is always at about 20 — 30 fs. 



criteria which can alter this trend since it will explicitly mod- 
ify the couplings. Then secondly, as noticed from the wave- 
function propagation, diffusion is faster along the x axis than 
it is along the y axis. This difference is largest for small values 
of disorder (about a 3-fold difference)) where the propagation 
is fastest, in the quasi ballistic regime. Going to larger disor- 
der, the wavefunction spreads much slower and propagation is 
reduced in both directions. 

Finally, we investigated the transport as a function of static 
disorder as well. The three dimensional surface plots for TC 
and TDBC are shown in Fig. 13. We notice that diffusion is 
strongly dependent on and much less dependent on a. For 
each fixed value of F"^, the largest variation of Da over the 
a interval is of about 30% of the largest values. On the other 
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Figure 11. Top panel: Density of exciton states as a function of en- 
ergy of states computed for an aggregate of 51x51 TDBC molecules 
and averaged over 1000 realizations of static disorder a — 70 meV. 
The zero of energy corresponds to the electronic transition of a sin- 
gle molecule. The energy of the central site (injection point) in the 
lattice is assumed to be zero. Bottom panel: Plot of second moments 
of the wave function at time t — 100 fs as a function of the initial 
injection energy for TDBC with static disorder ct = 70 meV and 
r''' = 30meV. 
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Figure 12. Comparison of diffusion constants as a function of 
dynamic disorder F'^ for TC, TDBC and U3 with static disorder 
a = 70meV. As dynamic disorder is increased, the diffusion con- 
stants decrease. Further, xx diffusion constants are larger than yy 
components for all values of dynamic disorder. Finally, diffusion is 
larger for U3 than it is for TDBC and TC. This can be explained by 
the molecular coupling parameters as described in the text. 



hand, the largest variation of Da at fixed a and varying 
is of the order of 80% of the largest value. Due to the pres- 
ence of static disorder the dependence of the diffusion coeffi- 
cient on the dephasing rate F"^ deviates from the conventional 
D (X l/F"^ derived within the Haken-Strobl-Reineker model 
for homogeneous systems [46]. 

Recently, Akselrod et. al. [22], have estimated the exciton 
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diffusion length in thin-film J-aggregates of TDBC molecules 
based on an exciton-exciton annihilation experiment at room 
temperature. In that study, the exciton lifetime obtained 
from time-dependent photoluminescence was determined to 
be Texp — 45 ps, and an expression for three-dimensional ex- 
citon diffusion was used to determine the annihilation rate. 
However, the sulphonated group side chains are about 6 A 
long and in addition, in the growth process a layer of poly- 
mer molecules is introduced between the J-aggregate layers. 
As a result the spacing between the monolayers of fluorescent 
dyes is several times larger than the distance between nearest 
neighbor molecules in a layer Therefore we can assume that 
diffusion in these aggregates is two dimensional rather than 
three dimensional. Using a 2D model, we find that the ex- 
perimental exciton diffusion length is about £cxp ~ 60 nm. 
We can estimate the exciton diffusion length along the i-th di- 
rection using the measured lifetime and a computed diffusion 
coefficient as 

= v/2A»Te,p. (23) 

For the value of the dynamic disorder F"^ =30 meV, which 
should approximately correspond to room temperature, as dis- 
cussed in Section 11 B, and the static disorder a = 70 meV we 
find £.j., w 200 nm and £y ~ 100 nm. These values are in a 
good qualitative agreement with the measured one. However, 
the quantitative discrepancy can be due to a number of factors 
such as different lattice constants of the aggregate, additional 
exciton dephasing and relaxation channels, and also domain 
boundaries in the experimental structures. All of these aspects 
can be incorporated into the model provided that one can ex- 
tract the actual parameters should from experiments. 

V. CONCLUSIONS 

In this article, a mixed model combining an open quan- 
tum systems approach to ab-initio calculations has been em- 
ployed to gain insight on the exciton dynamics of thin-film 
J-aggregates. This model can capture both coherent and in- 
coherent transport and allows for a detailed study of transport 
parameters such as diffusion coefficients and diffusion length. 
The role of the initial state of the system on transport can also 
be captured. Further one can investigate all these aspects as a 
function of the structure of the aggregates through the lattice 
parameters. 

As an example the model was applied to three different 
cyanine dye aggregates. Within this model we conclude that 
transport depends explicitly on the molecular properties of the 
monomers which compose the aggregate. In particular, for 
molecules packed in a brickstone arrangement, transport in- 
creases with the monomers' transition dipoles and hence with 
the coupling between monomers. Furthermore, the coupling 
induces a preferential direction for transport which leads to 
an anisotropic spread of populations and coherences. Such 
directionality is robust to both static and dynamic disorder 
within the investigated ranges and does not change for dif- 
ferent molecules. This model has permitted the identification 



of timescales for the different transport regimes. A ballistic 
regime is present for all values of disorder at times smaller 
than 20 fs while afterwards a diffusive regime is observed. 
The transport is also determined by the choice of the initial 
condition, at a specific value of injection energy a maximum 
in diffusion is observed. Investigation on the origin of the ex- 
act position of this maximum are planned but qualitatively it 
can be explained by the large DOS located between J-band 
and monomer transition. 

The obtained diffusion length is in good agreement with 
experimental results, however a more accurate comparison 
would be possible only including relaxation in the excited 
state manifold. The model does not account for exciton do- 
mains which would reduce diffusion, and further we cannot 
investigate the low temperature regime. Work in these direc- 
tions is currently in progress in our research group. 

The efficient exciton transport observed in these thin-film 
aggregates and the possibility of tuning this transport by 
the choice of monomers or by selecting the initial condition 
makes these aggregates good candidates for devices where 
large exciton diffusion lengths are sought. In particular cou- 
pling them to optical micro-cavities [96] opens the road to a 
range of control possibilities which we plan to study in future 
work. Coupling these large-exciton diffusion length materials 
to high-hole mobility materials [97] might also provide some 
advantages for technological appUcations such as all-organic 
photon detectors. 
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Appendix A: Single molecule Hamiltonian 

In the harmonic approximation, the single-molecule elec- 
tron excitation Hamiltonian can be written as 




(Al) 



where J7, is the energy of i-th electronic transition, b]^ and 
big are creation and annihilation operators for the g-th vi- 
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Figure 13. Diffusion coefficients of singlet excitons in 2D J-aggregates of TC (panel a) and TDBC (panel b) as functions of static o and 
dynamic disorders V" . The initial condition for propagation was injection at Ein — —281 meV for TDBC and Ei„ = — 287meV for TC. 
The red color corresponds to smaller diffusion while the yellow indicates larger diffusion coefficients. 



bration and is the frequency of the corresponding vibra- 
tional mode. Within the visible range of the spectrum, the 
lowest electronic transition in fluorescent dyes has the largest 
oscillator strength. This transition determines the formation 
of a J-band in molecular aggregates. An approximation in- 
volving only the two lowest electronic states - the ground \g) 
and the first excited |e) states of a monomer - can provide 
good description of the excitation dynamics as long asthe sec- 
ond excited state is far enough from the first state. This ap- 
proximation has been used previously to model excitons in 
J-aggregates [15, 98, 99]. In the two-level approximation Eq. 
Al simplifies to 



4„oi = fie |e) (e| + ^ Wqfejfe, + ^ |e) (e| (SJ + 
q q 

(A2) 

where we have assumed that the frequencies of the vibrational 
modes are the same for all electronic states, and the renor- 
malized electronic transition frequency fig — il^ + J2q \ 
involves the reorganization energies A, associated with each 
mode q. The exciton-vibration coupling coefficient is defined 
as Kg — Xq/Aq, where the boson operators associated with 
different electronic states are related by beq = bgq+Aq. In Eq. 
A2 all vibrational modes correspond to the electronic ground 
state and the corresponding ground state indices are omitted 
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Normal coordinate, Q 



Figure 14. Schematic energy profile of a single vibrational mode in a 
two-level molecule. A is a reorganization energy and AqIs a normal 
coordinate displacement. 



for simplicity. The phonon creation/annihilation operator dis- 
placement Aq is related to the normal coordinate displace- 
ment operator Ag^ as A^ = Ag^ • ^^uj^jh. A schematic 



diagram of a two-level molecule is shown in Fig. 14. 



Appendix B: Exciton manifolds 

For an aggregate consisting of iV molecules all 2^ — 1 ex- 
cited states can be classified into N different exciton mani- 
folds - single exciton X^^^ , bi-exciton X*^^', etc, see Fig. 15. 
The energy separation between the manifolds is of the same 
order of magnitude as that of the single exciton energy. The 
exiton manifolds are coupled by exciton relaxation or annihi- 
lation processes and also by interaction with external radia- 
tion. Within each subspace, the exciton states of monomers 
are coupled to each other by the Coulomb interaction, which 
gives rise to a resonant excitation transfer between molecules 
[43]. The same Coulomb interaction also renormalizes the 
energies of single monomer excitations. However, transitions 
between different exciton manifolds induced by the Coulomb 
interaction can be ignored. Thus, one can study single ex- 
citon dynamics by constraining the exciton space to the sin- 
gle exciton manifold X*^^' and the ground state X'^'^^ . In Fig. 
15 we show a schematic diagram of exciton manifolds to- 
gether with the allowed transitions for a system of three cou- 
pled molecules. Unlike molecular crystals, delocalized exci- 
ton bands are not formed in molecular aggregates due to a 
strong competition between the Coulomb interaction, which 
tends to delocalize the exciton and to the exciton-vibration 
coupling and structure disorder, which favor localization. 
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